Introduce types, e.g. `DivDiffData`, to get around complicated broadcast machinery #142

DanielVandH · 2025-10-21T10:36:05Z

@dlfivefifty Is this roughly what you had in mind for trying to get around some of the broadcast machinery? I'm simply caching the broadcast results in a CachedBroadcastVector with cache_filldata! implemented as

function LazyArrays.cache_filldata!(K::CachedBroadcastVector{T, F, A, B}, inds) where {T, F, A, B}
    @show inds, K.datasize
    f, a, b = K.f, K.A, K.B
    resizedata!(a, maximum(inds)) 
    resizedata!(b, maximum(inds))
    @inbounds for k in inds
        K.data[k] = f(a[k], b[k])
    end
end

It's currently only implemented in divdiff to get an idea of what it looks like:

function divdiff(Q::SemiclassicalJacobi, P::SemiclassicalJacobi)
    A,B,_ = recurrencecoefficients(P)
    α,β,_ = recurrencecoefficients(Q)

    d = AccumulateAbstractVector(*, CBV(/, A, Vcat(1,α)))
    v1 = AccumulateAbstractVector(+, CBV(/, B, A))
    v2 = MulAddAccumulate(CBV(/, Vcat(0,0,α[2:∞]), α), CBV(/, Vcat(0,CBV(/, β,  α)),  α))
    v3 = AccumulateAbstractVector(*, Vcat(A[1]A[2], CBV(/, A[3:∞], α)))
    op1 = CBV(*, 1:∞, d)
    op2 = CBV(/, B[2:end], A[2:end])
    op3 = CBV(+, v1, op2)
    op4 = CBV(*, 1:∞, op3)
    op5 = CBV(/, β, α)
    op6 = CBV(*, α, v2)
    op7 = CBV(+, op5, op6)
    op8 = CBV(*, 2:∞, op7)
    op9 = CBV(-, op4, op8)
    op10 = CBV(*, op9, v3) 
    return _BandedMatrix(Vcat(op1', op10'), ∞, 2,-1)'
    # return _BandedMatrix(Vcat(((1:∞) .* d)', (((1:∞) .* (v1 .+ B[2:end]./A[2:end]) .- (2:∞) .* (α .* v2 .+ β ./ α)) .* v3)'), ∞, 2,-1)'
end

which is a bit annoying as it requires a bunch of individual CBV calls. I did initially start with a CachedFusedBroadcastVector but it's somewhat unreadable as I need to manually build the Tuple of functions and vectors. I also tried specific structs for each type of operation, e.g. RangeTimesAccumulate, but getting this type of struct correct seems nicer.

There's an issue with this implementation though that I'm not sure I understand. It keeps trying to indefinitely resize the vector in some cases. For example, putting a @show inds, K.datasize inside LazyArrays.cache_filldata!(K::CachedBroadcastVector{T, F, A, B}, inds) where {T, F, A, B} and differentiating a weighted basis,

julia> diff(Weighted(SemiclassicalJacobi(2.0, 2, 2, 2)))
...
[truncated]
...
(inds, K.datasize) = (454:454, (453,))
(inds, K.datasize) = (455:455, (454,))
(inds, K.datasize) = ERROR: InterruptException:
Stacktrace:
  [1] try_yieldto(undo::typeof(Base.ensure_rescheduled))
    @ Base .\task.jl:1128
  [2] wait()
    @ Base .\task.jl:1200
  [3] uv_write(s::Base.TTY, p::Ptr{UInt8}, n::UInt64)
    @ Base .\stream.jl:1081
  [4] unsafe_write(s::Base.TTY, p::Ptr{UInt8}, n::UInt64)
    @ Base .\stream.jl:1154
  [5] write
    @ .\strings\io.jl:246 [inlined]
  [6] print
    @ .\strings\io.jl:248 [inlined]
  [7] print(::Base.TTY, ::String, ::String, ::String)
    @ Base .\strings\io.jl:46
  [8] println(x1::String, x2::String)
    @ Base .\coreio.jl:6
  [9] cache_filldata!(K::SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{…}, inds::UnitRange{…})
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\derivatives.jl:52
 [10] _vec_resizedata!(B::SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{…}, n::Int64)
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\cache.jl:246
 [11] resizedata!
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\cache.jl:254 [inlined]
 [12] resizedata!
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\cache.jl:222 [inlined]
 [13] getindex
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\cache.jl:79 [inlined]
 [14] getindex
    @ C:\Users\djv23\.julia\juliaup\julia-1.12.1+0.x64.w64.mingw32\share\julia\stdlib\v1.12\LinearAlgebra\src\adjtrans.jl:342 [inlined]
 [15] copyto_unaliased!(deststyle::IndexLinear, dest::LinearAlgebra.Adjoint{…}, srcstyle::IndexLinear, src::LinearAlgebra.Adjoint{…})
    @ Base .\abstractarray.jl:1090
 [16] copyto!(dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ Base .\abstractarray.jl:1070
 [17] copyto!_layout(::ArrayLayouts.DualLayout{…}, ::LazyArrays.LazyLayout, dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ArrayLayouts.jl:266
 [18] copyto!_layout(dlay::ArrayLayouts.DualLayout{…}, ::ArrayLayouts.DualLayout{…}, dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\memorylayout.jl:310
 [19] copyto!_layout(dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ArrayLayouts.jl:272
 [20] copyto!(dest::LinearAlgebra.Adjoint{…}, src::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ArrayLayouts.jl:285
 [21] copymutable
    @ .\abstractarray.jl:1201 [inlined]
 [22] copy
    @ .\abstractarray.jl:1144 [inlined]
 [23] map
    @ .\tuple.jl:359 [inlined]
 [24] copy(f::ApplyArray{Float64, 2, typeof(vcat), Tuple{LinearAlgebra.Adjoint{…}, LinearAlgebra.Adjoint{…}}})
    @ LazyArrays C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\lazyconcat.jl:97
 [25] copy(B::BandedMatrices.BandedMatrix{Float64, ApplyArray{…}, InfiniteArrays.OneToInf{…}})
    @ BandedMatrices C:\Users\djv23\.julia\packages\BandedMatrices\Vz6Ql\src\banded\BandedMatrix.jl:167
 [26] copy(A::LinearAlgebra.Adjoint{Float64, BandedMatrices.BandedMatrix{Float64, ApplyArray{…}, InfiniteArrays.OneToInf{…}}})
    @ InfiniteArraysBandedMatricesExt C:\Users\djv23\.julia\packages\InfiniteArrays\LdDTO\ext\InfiniteArraysBandedMatricesExt.jl:416
 [27] _copy_oftype(A::LinearAlgebra.Adjoint{Float64, BandedMatrices.BandedMatrix{…}}, ::Type{Float64})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ArrayLayouts.jl:70
 [28] copy(M::Mul{…})
    @ LazyArraysBandedMatricesExt C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\ext\LazyArraysBandedMatricesExt.jl:571
 [29] materialize
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\mul.jl:137 [inlined]
 [30] mul
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\mul.jl:138 [inlined]
 [31] *(A::LinearAlgebra.Diagonal{…}, B::LinearAlgebra.Adjoint{…})
    @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\mul.jl:303
 [32] _copy_ldiv_mul
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\linalg\inv.jl:115 [inlined]
 [33] copy
    @ C:\Users\djv23\.julia\packages\LazyArrays\gHXdm\src\linalg\inv.jl:119 [inlined]
 [34] copy
    @ C:\Users\djv23\.julia\packages\ContinuumArrays\TDola\src\bases\bases.jl:368 [inlined]
 [35] copy(L::Ldiv{…})
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\SemiclassicalOrthogonalPolynomials.jl:416
 [36] materialize
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ldiv.jl:22 [inlined]
 [37] ldiv
    @ C:\Users\djv23\.julia\packages\ArrayLayouts\FxFFC\src\ldiv.jl:98 [inlined]
 [38] \(A::SemiclassicalJacobi{Float64}, B::QuasiArrays.ApplyQuasiMatrix{Float64, typeof(*), Tuple{…}})
    @ QuasiArrays C:\Users\djv23\.julia\packages\QuasiArrays\wlEWh\src\matmul.jl:34
 [39] divdiff(wP::Weighted{Float64, SemiclassicalJacobi{Float64}}, wQ::Weighted{Float64, SemiclassicalJacobi{Float64}})
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\derivatives.jl:110
 [40] diff(wQ::Weighted{Float64, SemiclassicalJacobi{Float64}}; dims::Int64)
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\derivatives.jl:115
 [41] diff(wQ::Weighted{Float64, SemiclassicalJacobi{Float64}})
    @ SemiclassicalOrthogonalPolynomials c:\Users\djv23\.julia\dev\SemiclassicalOrthogonalPolynomials.jl\src\derivatives.jl:113
 [42] top-level scope
    @ REPL[19]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> diff(Weighted(SemiclassicalJacobi(2.0, 2, 2, 2)))

Any idea what's going on with this implementation that it keeps wanting to resize? I thought _vec_resizedata!(B::AbstractVector, n) should get around this automatically. Note that the indexing does work sometimes though, e.g.

julia> diff(SemiclassicalJacobi(2.0, 2, 2, 2)).args[2]
ℵ₀×ℵ₀ adjoint(::BandedMatrices.BandedMatrix{Float64, ApplyArray{Float64, 2, typeof(vcat), Tuple{LinearAlgebra.Adjoint{Float64, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), InfiniteArrays.InfUnitRange{Int64}, Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}, LinearAlgebra.Adjoint{Float64, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(-), SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), InfiniteArrays.InfUnitRange{Int64}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(+), Cumsum{Float64, 1, AbstractVector{Float64}}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(/), BroadcastVector{Float64, typeof(-), Tuple{BroadcastVector{Float64, typeof(/), Tuple{SubArray{Float64, 1, InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}, SubArray{Float64, 1, InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}}}, BroadcastVector{Float64, typeof(inv), Tuple{SubArray{Float64, 1, InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}, Tuple{InfiniteArrays.InfUnitRange{Int64}}, false}}}}}}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), InfiniteArrays.InfUnitRange{Int64}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(+), SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(/), BroadcastVector{Float64, typeof(-), Tuple{BroadcastVector{Float64, typeof(/), Tuple{InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}, InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}}}}}, BroadcastVector{Float64, typeof(inv), Tuple{InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}}}}, SemiclassicalOrthogonalPolynomials.CachedBroadcastVector{Float64, typeof(*), BroadcastVector{Float64, typeof(inv), Tuple{InfiniteLinearAlgebra.TridiagonalConjugationBand{Float64}}}, MulAddAccumulate{Float64}}}}}, Accumulate{Float64, 1, typeof(*), Vector{Float64}, AbstractVector{Float64}}}}}}, InfiniteArrays.OneToInf{Int64}}) with eltype Float64 with indices OneToInf()×OneToInf():
  ⋅   5.40256  -1.09972    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅       …  
  ⋅    ⋅        8.0885   -1.72961    ⋅         ⋅         ⋅         ⋅         ⋅
  ⋅    ⋅         ⋅       10.4884   -2.26406    ⋅         ⋅         ⋅         ⋅
  ⋅    ⋅         ⋅         ⋅       12.7791   -2.74673    ⋅         ⋅         ⋅
  ⋅    ⋅         ⋅         ⋅         ⋅       15.0146   -3.19782    ⋅         ⋅
  ⋅    ⋅         ⋅         ⋅         ⋅         ⋅       17.2173   -3.62826    ⋅       …
  ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅       19.3988   -4.04445
  ⋅    ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅       21.5657
 ⋮                                            ⋮                                      ⋱

works and uses CachedBroadcastVector.

dlfivefifty · 2025-10-21T13:08:15Z

What's the difference between a CachedBroadcastVector and a CachedVector{<:Any,<:BroadcastVector}?

We definitely do not want to cache every single operation

DanielVandH · 2025-10-21T14:24:45Z

I think the issue was with type inference for nested BroadcastVectors (or maybe just for caches of them) gives an eltype of Union{}.

Maybe if instead of using a cache for all the individual broadcasts I just build up all the BroadcastVectors and cache those

dlfivefifty · 2025-10-21T15:39:59Z

No the issue has nothing to do with type inference. It' that if you broadcast cached vectors then call to getindex have no way to tell the underlying vectors to call resizedata!. And we don't have the ability to reduce a getindex to a broadcast of the underlying cached data.

The proposal was to do a single type to wrap and cache the entire vector, not to introduce more generic lazy arrays that are composed.

DanielVandH · 2025-10-21T17:58:34Z

I think this is closer to what you want? I still need to figure out why e.g. (taken directly from the tests)

t = 2
P = SemiclassicalJacobi(t, -0.5, -0.5, -0.5)
Q = SemiclassicalJacobi(t,  0.5,  0.5,  0.5, P)
x = axes(P,1)
D = Derivative(x)
D = Q \ (D*P)

tries to indefinitely resize the cache though.

Rough layout

afdba25

Remove CBV

ae38024

DanielVandH changed the title ~~Implement CachedBroadcastVector to get around broadcast machinery~~ Introduce types, e.g. DivDiffData, to get around complicated broadcast machinery Oct 21, 2025

Add more datas for the HalfWeighted derivatives

c5c00b1

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Introduce types, e.g. `DivDiffData`, to get around complicated broadcast machinery #142

Introduce types, e.g. `DivDiffData`, to get around complicated broadcast machinery #142

DanielVandH commented Oct 21, 2025

Uh oh!

dlfivefifty commented Oct 21, 2025

Uh oh!

DanielVandH commented Oct 21, 2025

Uh oh!

dlfivefifty commented Oct 21, 2025

Uh oh!

DanielVandH commented Oct 21, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Introduce types, e.g. DivDiffData, to get around complicated broadcast machinery #142

Are you sure you want to change the base?

Introduce types, e.g. DivDiffData, to get around complicated broadcast machinery #142

Conversation

DanielVandH commented Oct 21, 2025

Uh oh!

dlfivefifty commented Oct 21, 2025

Uh oh!

DanielVandH commented Oct 21, 2025

Uh oh!

dlfivefifty commented Oct 21, 2025

Uh oh!

DanielVandH commented Oct 21, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Introduce types, e.g. `DivDiffData`, to get around complicated broadcast machinery #142

Introduce types, e.g. `DivDiffData`, to get around complicated broadcast machinery #142